Instability of the chiral d-wave RVB state for the Heisenberg model on triangular 
lattice and an improvement of the Gutzwiller approximation 
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Through Variational Monte Carlo simulation we show the d-wave RVB pairing in the Heisenberg 
model on triangular lattice can be better described in terms of a two component order parameter. 
The fully gapped chiral d-wave RVB state, which is predicted by the mean field theory to be the 
unique minimum of variational energy in the two dimensional representation space of d-wave pairing, 
is found to be actually a local maximum and the true minimum of energy is reached by the non- 
chiral d xy state with line nodes. We also find that the usual Gutzwiller approximation, which enjoys 
great success for the square lattice system, fails badly on the triangular lattice as a result of the 
geometric frustration inherent of the system. An improved version of the Gutzwiiler approximation 
is proposed and is found to give a much better results than the usual one. 

PACS numbers: 



The search for spin liquid ground state on geometri- 
cally frustrated quantum magnet is a central issue in the 
study of the strongly correlated electron systems. On the 
one hand, the spin liquid state represents a novel state of 
matter beyond the Landau-Ginzburg description. On the 
other hand, the study of the spin liquid state is closely re- 
lated to the study of exotic superconductivity in strongly 
correlated electron system. The most famous example in 
this respect is the high-T c cuprates. Although the parent 
compounds of the cuprates(an ideal quantum antiferro- 
magnet described by the Heisenberg model on square lat- 
tice) exhibit magnetic order, the study of spin liquid state 
for the Heisenberg model on square lattice nevertheless 
predicts unambiguously the d-wave paring symmetry for 
the superconducting order at finite doping and help to 
clarify the nature of the pseudogap phase l]. 

The discovery of superconductivity in hydrated 
Na x Co02[£l triggered another round of intensive re- 
searches in this field as a result of the geometric frus- 
tration and the strong electron correlation inherent of 
the system which has a triangular lattice [3j-|6|]. Histor- 
ically, it is first on the triangular lattice that the very 
concept of spin liquid state was first proposed [7]. How- 
ever, after many year's of intensive investigation the na- 
ture of the superconducting state remains illusive. While 
NMR measurement suggests a spin singlet pairing with 
line node, theories of the t-J model on triangular lattice 
predicts a chiral d-wave state with full gap and time re- 
versal symmetry breaking. In recent years, the possibility 
of spin liquid state and exotic superconducting state(with 
self-doped charge carrier) in half-filled systems with an 
anisotropic triangular lattice and reduced strength of lo- 
cal correlation also arose great interest in the field [8]. 
Here the same uncertainty on the nature of the spin liq- 
uid state exist |9l. llCll. Il2i? ]. While some theory predicts 
a spin liquid state with a open Fermi surface for spin ex- 
citation, others argue spin pairing is unavoidable in the 
spin liquid state. 



As compared to the square lattice, the triangular lat- 
tice is more complex in that there are more choices for the 
spin pairing pattern on the triangular lattice. According 
to group theory, spin pairing between neighboring sites 
on the square lattice can be either extended s-wave or 
d-wave, both of which belong to a one dimensional ir- 
reducible representation of the C^ v group. However, on 
the triangular lattice with a point group of Cq v , the d- 
wave pairing channel becomes two dimensional and it is a 
non-trial problem to decide what kind of pairing pattern 
is actually realized within the two dimensional represen- 
tation space of the d-wave pairing(the extended s-wave 
pairing is in general less stable than d-wave pairing for 
system with strong local correlation). 

The slave Boson mean field theory and the Gutzwiller 
approximation based on it is widely adopted to study the 
problem of spin liquid and the pairing symmetry of the 
strongly correlated system [1]. Such a approach receives 
great succuss on the square lattice. For example, it pre- 
dicts correctly the d-wave pairing symmetry of the high- 
Tc cuprates 13]. With the Gutzwiller approximation, the 
mean field theory can even give a quantitative correct es- 
timation of the variational energy with an error less than 
5 percent [lj, Ilia]- On the triangular lattice, the slave bo- 
son mean field theory predicts unambiguously that the 
chiral d-wave pairing with time reversal symmetry break- 
ing and a full gap is the most stable in two dimensional 
representation space and its condensation energy is more 
than 30 percent lower than that of the non-chiral(real) 
d-wave pairing state. The chiral d-wave state is widely 
believed to be the best choice for the triangular lattice. 

In this paper, we show the usual Gutzwiller approxi- 
mation fails badly for the Heisenberg model on the tri- 
angular lattice. It overestimates the condensation energy 
by more than 50 percent and predicts incorrect order for 
the relative stability of the chiral d-wave and non-chiral 
d-wave state. Through direct variational Monte Carlo 
simulation, we find the chiral d-wave RVB state is actu- 



ally a local maximum in the two dimensional represen- 
tation space of d-wave pairing and the true minimum of 
variational energy is reached by the non-chiral d xy state 
with line nodes. The anisotropy of the variational energy 
in the representation space is found to be much smaller 
than that predicted by the mean field theory and the 
RVB pairing should be more appropriately described in 
terms of a two component order parameter. We also find 
the failure of the usual Gutzwiller approximation on tri- 
angular lattice can be attributed to the geometric frus- 
tration inherent of the system. Based on this understand- 
ing, we propose an improved Gutzwiller approximation 
scheme and find it works much better than the usual one. 
The Hciscnbcrg model under consideration reads 



H 



j 2_^ Si • Sj 



(i) 



<hj> 



here X}<j 7 > means sum over nearest neighbors on the tri- 
angular lattice. The slave Boson mean field theory of the 
RVB state is established by introducing the slave parti- 
cle representation of the spin operator Si = 5// a &a,pfi,p 
and then use the mean field order parameter Xij = 

(ftafj-a) and A y = {za,pfi t0l f] t p) to decouple the Hamil- 
tonian written in terms of the slave particles. The slave 
particles must be subjected to the no double occupancy 
constraint to represent the spin operator faithfully. In 
mean field treatment, such local constraint is relaxed to 
a constraint on the average particle number on each site. 
Assuming translational invariance for the RVB order 
parameter \i,j an( i A,,j, the mean field Hamiltonian for 
the slave particle has the form, 



HmF — y^£kf k , a fk,a + /(Afc/, 



k . a 

whose ground state reads 



t A 



\BCS)=l[(u k + v k fyi H )\0), 



(3) 



in which — 



Ek = \/£,l + A£ . Instead of deter- 
i_2 from the self-consistent equation, 



mining \i,j an d A 

we will take the wave function Eq.(3) as a variational 
description of the system at the mean field level. For 
simplicity, Xi,j wm be assumed to be a real constant and 
be nonzero for nearest neighboring sites only. The hop- 
ping energy £j. then reads 



j. c\ / i n X V 'J'^U \ 

£k = — 2(cos k x + 2 cos — cos - 1 



//, 



(4) 



here we have set the hopping integral as one. The chem- 
ical potential is determined by the mean field equation 
for the particle number and is not treated as an indepen- 
dent variational parameter. Thus the magnitude of the 
gap function, A, is the only variational parameter in the 
theory. 



The physical wave function of the spin system should 
satisfy the no double occupancy. This can be achieved by 
taking the Gutzwiller projection on the mean field wave 
function. The Gutzwiller projected mean field state takes 
the form 



I*) = Pa 
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(5) 



where the Cooper pair wave function a,i-j is given by 

A k 
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ik(Ri-Rj 



\ (6) 



and Pq is the projection into the subspace of no dou- 
ble occupancy. The Gutzwiiler projected wave function 
is usually treated by the variational Monte Carlo simula- 
tion method as analytical calculation is difficult for it. An 
estimate of variational energy in the Gutzwiller projected 
wave function can be given by the so called Gutzwiller ap- 
proximation, in which the expectation value of a physical 
observable in the Gutzwiller projected state is approxi- 
mated by its expectation value in the mean field state 
multiplied with a Gutzwiller factor. For the Heisenberg 
system, the Gutzwiller factor for spin correlation is usu- 
ally taken to be g s = 4. 

Now we determine the structure of the pairing order 
parameter. Here we will consider the d-wave pairing only. 
On the triangular lattice, the d-wave pairing belongs to 
a two dimensional irreducible representation of the point 
group of Cq v and is thus more complex than its coun- 
terpart on the square lattice, which has a unique basis 
function of the form of d x 2_ y 2. The two basis functions 
(2) of this two dimensional representation have the form of 
d x 2_ y 2 and d xy . The symmetry of the system guarantee 
the invariance of the free energy under rotation in the 
two dimensional representation space at the quadratic 
level. Any anisotropy in the representation space must 
be attributed to higher order terms in the Landau expan- 
sion. Thus, if such anisotropy is small, one should better 
describe the pairing state with a multi-component order 
parameter. 

In momentum space, these two basis functions of 
the d-wave representation have the form of $i(k) = 
cos k x - cos \ cos ^ke. an d $ 2 (k) = V3 sin \ sin ^|^. 
In real space, they behave as cos2iy9 and sm2ip respec- 
tively(see Fig.l), where ip is the angle the bond made 
with a;— axis. The general order parameter is given by 
Ak = A[r/i$i(k) + ?72 < i ) 2(k)], where 771 and 772 are two 
complex numbers satisfying |7yi| 2 + \r)2\ 2 — 1. Up to a 
global phase, the gap function can then be parameter- 
ized as follows 

A k = A[cos0$i(k) +sin0exp(i0)$ 2 (k))], (7) 

in which A, 9 and <f> are all real numbers. In this param- 
cterizition, the chiral d-wave state is given by {9,<j>) — 
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the chiral d x 2_ y 2 +id xy state and the non-chiral d xy state 
as a function of A. It is found that the chiral d-wave state 
enjoys a more than 30 percent lower condensation energy 
than its non-chiral counterpart. For this reason, the two 
component nature of the order parameter is irrelevant 
and we can take the chiral d-wave state as the unique 
choice for spin pairing. 



FIG. 1: The two basis functions of d-wave pairing on trian- 
gular lattice. (a,)d x 2_ y 2,(b)d xy . 



(■j, ?). Owing to the gauge symmetry and symmetry of 
space inversion and time reversal, it is easy to see that we 
only need to consider the value of 9 and <f> in the region 

(M)e[0,7r/2]x[0,7r/2]. 

At the mean field level, the chiral d x 2_ y 2 + id xy state 
gives the lowest energy in this space while the maximum 
in energy is reached by the non-chiral d xy state. Note 
also that the energy is periodic in 9 when cj> — with a 
period of 7r/3, as a result of the six- fold rotational sym- 
metry of the triangular lattice. A full scan of the vari- 
ational energy(optimized with respect to the parameter 
A) estimated from the usual Gutzwiller approximation 
with g s = 4 is shown in Fig. 2. For any given 9, a real 
gap function always gives higher variational energy than 
a complex gap function. 
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FIG. 2: A full scan of the variational energy (optimized with 
respect to A)calculated from the Gutzwiller approximation 
with g s — 4 in the representation space of d-wave pairing. 

Figure 3 provides a more quantitative measure on the 
anisotropy in the representation space at the mean field 
level, in which the variational energy calculated from the 
Gutzwiller approximation with g s = 4 is compared for 
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FIG. 3: The variational energy as a function of A calculated 
from the Gutzwiller approximation with g s = 4 for the chiral 
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u 2 + id xy state and the non-chiral d xy state. 



Now we present the results of Variational Monte Carlo 
simulation. The simulation is done on a 14 x 14 lattice 
with periodic-antiperiodic boundary condition. 10 sta- 
tistically independent samples are used to estimate the 
variational energy. The result for the chiral d x 2_ y 2 +id xy 
state and the non-chiral d xy state are shown in Figure 4. 
The result differs from the mean field result in three im- 
portant aspects. Firstly, the order of relative stability 
between the chiral d x 2_ y 2 +id xy state and the non-chiral 
d xy state is now reversed. Secondly, the difference in con- 
densation energy is much smaller than that predicted by 
the mean field theory, being below 3 percent. Lastly, the 
variational energy estimated from the Gutzwiller approx- 
imation is about 50 percent lower than that calculated 
directly from the Monte Carlo simulation. We thus con- 
clude that the usual Gutzwiller approximation fails badly 
on the triangular lattice. 

A full scan of the variational energy(optimized with re- 
spect to A) as a function of 9 and <j> is shown in Figure 5. 
The simulation is done on a 16 x 16 lattice with periodic- 
antiperiodic boundary condition. The variational energy 
is approximately periodic in 9 with a period of 7r/3 when 
4> = 0, as required by the six- fold rotational symmetry of 
the triangular lattice. The small deviation from such 
a periodicity, especially the small extra peaks around 
(9,<j>) = (7r/12,0) and {6,4>) = (5tt/12,0), are caused 
by finite size effect which is the strongest for the non- 
chiral state with line node. From the Figure we see the 
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FIG. 4: The variational energy as a function of A for the chiral 
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+ id xy state and the non-chiral d xy state calculated 



from VMC simulation. The simulation is done on a 14 x 14 
lattice with periodic-antiperiodic boundary condition. 



FIG. 6: The variational energy of the chiral d x i_ y i + id xy 
state(thick lines) and the non-chiral d xy state(thin lines) 
as a function of A calculated from VMC simulation (black 
lines), one site Gutzwiller approximation(red lines), two site 
Gutzwiller approximation (green lines) and the three site 
Gutzwiller approximation (blue lines). 



state as follows, 




FIG. 5: A full scan of the variational energy (optimized with 
respect to A) in the representation space of d-wave pairing 
calculated from VMC simulation. The simulation is done on a 
16 x 16 lattice with periodic-antiperiodic boundary condition. 



true minimum of the variational energy is reached by the 



non-chiral d xy state, while the chiral d 



now becomes a local maximum. At the same time, the 
anisotropy of the condensation energy in the represen- 
tation space, being less than 5 percent, is much smaller 
than that predicted by the mean field theory. For this 
reason, it is more appropriate to describe to the d-wave 
pairing in terms of a two component order parameter. 

Then, why does the usual Gutzwiller approximation, 
which perform well for the square lattice, fails so badly 
on the triangular lattice? To understand this, we rewrite 
the expectation value (sfsj) in the Gutzwiller projected 
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(a)s z Ja), 



(8) 
in which a denotes one of the four spin configurations^! 
, H, 11, ||) on site i and j, W a and W' denote the prob- 
abilities for the spin configuration a to appear in the 
projected and mean field state. If we neglect the a de- 



pendence of the ratio tK$ , we get 



W> 



in which g s 



.E^(«M(°) 



QsiSiS^o, 



(9) 



L tot a 

in the mean field state 



()o denotes the expectation value 



We note that such a approx- 
is indeed a inde- 
pendent. In the BCS mean field state studied in this 
paper, the Gutzwiller factor g is easily calculated to be 



imation scheme becomes exact if |Kj 



?iven by g s = — 



-, where Xij = E CT 



C i,a c j,v)0, 



1 ^o- This is nothing but the two 



site version of the(improved) Gutzwiiler approximation 
derived by Hsu|15|. The usual Gutzwiller approximation 
with g s = 4 can be derived by simply neglecting the inter- 
site correlation and will be called the single site version 
of Gutzwiller approximation in the following. In Figure 
6, we show the variational energy calculated from both 
the two-site and the single site version of Gutzwiiler ap- 
proximation and compared them to the result of VMC 
simulation. Although the two-site Gutzwiller approxi- 
mation improve the result of the single site Gutzwiller 
approximation, none of them is truly satisfactory. 

The failure of the above approximation schemes can 
be attributed to the assumption that the ratio ||§ being 



a independent. On the triangular lattice, a given pair 
of nearest neighboring sites i and j are neighbored by a 
third site k. In the mean field state, the ratio between the 

probabilities for parallel and antiparallel spin alignment 

WS, 
on site i and j, namely R = ttttt, depends crucially on 

V TT 

the occupation of their common neighbor of site k. For 
example, if site k is empty, then site i and j can take the 
full advantage of forming singlet pair without been frus- 
trated by the spin on site k. Similarly, if site k is doubly 
occupied, then two antiparallel spins on site i and j can 
each form singlet pair with the electron on site k. On the 
other hand, if site k is singly occupied, no matter what 
the spin it has, its coupling to the two antiparallel spins 
on site i and j is frustrated. For this reason, the ratio R 
is larger when site k is cither empty or doubly occupied 
than it is singly occupied. This explains the breakdown 
of the assumption that J& being a independent. To sub- 
stantiate these arguments, we plot in Figure 7 the ratio 
R for the chiral d x 2_ y 2 + id xy state as a function of A 
when site k is either empty(i? e ), doubly occupied(i?d), 
or singly occupied(i? s ). From the figure we see R e is 
always higher than R s and Rd grows much faster than 
R s and exceeds it at large A. All these observations are 
consistent with the qualitative arguments raised above. 

The above reasoning also suggest a way to improve the 
Gutzwiller approximation. To reach this goal, we simply 
extend the two site version of the approximation to a 
three site version and use a to denotes the eight possible 
spin configurations on a elementary triangle of the lattice. 
In this scheme, the variational energy per triangle reads 



<(« 



-sfsl + 



s-s z k )) = g s {{s 



'^ + «M + *?aj0>o, (10) 



in which g s 



, j, k denotes the three neighbor- 



ing sites on an elementary triangle. The calculation of 
W° is straightforward but tedious and we will not present 
analytical expression for it here. For illustrative pur- 
pose, we show here the calculation of W® for one of the 
eight spin configuration, namely the configuration with 
all three spins aligned up. It is given by 



W° n = ((1 



n^)n^(l - Uj4.)rijt(l - n ki )n k ^) , (11) 



which can be evaluated with the Wick expansion. The 
expression for other W® can be similarly constructed and 
calculated. 

In Figure 6, we compare the variational energy calcu- 
lated from the one site, two site and the three site version 
of the Gutzwiiler approximation with the VMC results. 
While the two site version of Gutzwiller approximation 
is seen to have only a very limited improvement over the 
one site approximation, the three site approximation pro- 
vides a much more substantial improvement over the two. 
We thus conclude that the geometric frustration is indeed 
at the root of the failure of the usual Gutzwiller approx- 
imation on the triangular lattice system. It should also 
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FIG. 7: The ratio between the probabilities for antiparallel 
and parallel spin alignment on a pair of nearest neighboring 
sites i and j on the triangular lattice when their common 
neighbor site k is doubly occupied(-Rd), empty(7? e ) and singly 
occupied(i? s ). The calculation is done for the chiral d x 2_ y 2 + 
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be noted that the order of relative stability between the 
chiral and the non-chiral d-wave state is still incorrect 
at the three site version of the Gutzwiller approximation 
and higher order of approximation is needed to recover 
the correct order. 

In conclusion, we have studied the problem of RVB 
pairing symmetry for the Heisenberg model on the tri- 
angular lattice, which is relevant for both the study of 
nature of the Mott insulating state of the half filled Hub- 
bard model on the triangular lattice and the supercon- 
ducting state in the doped system. Unlike the square 
lattice system, the d-wave pairing on the triangular lat- 
tice require a description in terms of a two component 
order parameter. Contrary to the previous study based 
on usual Gutzwiiler approximation, we find the chiral- 
d wave state is actually a local maximum rather than 
local minimum in the two dimensional representation 
space of the d-wave pairing. We also find the anisotropy 
of condensation energy in the representation space is 
very small(less than 5 percent) and the true minimum 
is reached by the non-chiral d xy with line nodes. We 
find the usual Gutzwiller approximation, which perform 
well on square lattice, fails badly on the triangular lat- 
tice. The failure of the usual Gutzwiller approximation 
scheme is traced back to the geometric frustration inher- 
ent of the triangular lattice and an improved version of 
the Gutzwiller approximation is proposed based on this 
understanding and is found to result in substantial im- 
provement over the usual approximation scheme. 
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